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Abstract 

We study a two-dimensional Bose-Hubbard model at zero temperature with random local po- 
tentials in the presence of either uniform or binary disorder. Many low-energy metastable config- 
urations are found with virtually the same energy as the ground state's. These are characterized 
by the same blotchy pattern of the — in principle complex — non-zero local order parameter as the 
ground state. Yet, unlike the ground state, each island exhibits an overall random independent 
phase. The different phases in different coherent islands could provide a further explanation for 
the lack of coherence observed in experiments on Bose glasses. 

PACS numbers: 03.75.Lm, 05.30.Jp, 64.60.Cn 
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I. INTRODUCTION 



The problem of interacting bosons hopping across the sites of a disordered lattice has 
been addressed since the seminal paper on the Bose-Hubbard (BH) model by Fisher and 
co-workers There it was observed that the presence of disorder (in the form of random 
on-site potentials) enriches the phase diagram of the homogeneous model adding one further 
Bose-glass (BG) phase to the superfluid (SF) and Mott-insulator (MI) phases arising from tjhe 
competition between on-site repulsion and kinetic energy. Several studies followed Ref. 



11, 



employing various techniques and possibly addressing different realizations of disorder 



latter include models with diagonal disorder in one- dimensional (ID) lattices 



7|, |8| , and in 2D lattices |9|] , 
in 2D lattices Q, 



10l |. as well as off-diagonal disorder in ID lattices 




Recently, the interest in this subject has been boosted by the impressive advances in the 
field of cold atom trapping. Indeed, the BH model can be realized in terms of ultracold 
bosonic atoms trapped in optical lattices where disorder can be engineered with different 
experimental techniques. A speckle field can be superimposed otherwise an otherwise regu- 



lar lattice, resulting in random modulations of local lattice parameters [l^, [l^. A similar 
(pseudo)randomness can be obtained by superimposing two optical lattices with uncom- 
mensurate lattice constants Also, a regular optical lattice can be loaded with a 

small amount of a second (possibly fermionic) species with reduced mobility resulting from 
a quench in their hopping amplitude. Due to their interaction with the bosonic species, the 



atoms of the second species act as randomly localized scatterers 



It has been shown that a mean-field Gutzwiller a ppr oach captures the phases of the 
Bose-Hubbard model also in the presence of disorder [21|, |22| 



23l . l24j . Here we show that 



such an approach also highlights some features of the BG which support the reference to 
glasses in the name of such a phase. Indeed, the phase diagram of the disordered BH model 
is the subject of many investigations where the BG is usually characterized in terms of 
gaplessness, compressibility, superfluid and condensate fraction. Conversely, a few works 
focus on the typical features characterizing a glassy system 25|], such as, in particular, the 
exixtence of a complex energy landscape with a large number of local minima which should 
play a fundamental role in the relaxation dynamics. We expose the existence of many low- 



lying metastable states which differ from the BG ground-state essentially only for the phase 
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pattern of the order parameter. More specifically the BG phase is characterized by patches 
of finite local order parameter separated by regions where the local order parameter vanishes. 



This speckled pattern of the 
a vanishing superfluidity 



ocal order parameter results in a finite condensate fraction and 
23|. As we are going to illustrate, many states exist having the 
same distribution of local densities and (absolute value of) order parameters. The lowest 
energy is attained by the configuration where the phases of all the local order parameters 
are aligned along the same (arbitrary) direction. However, it turns out that many extremely 
stable configurations exist where the phase is aligned within each order-parameter patch, 
but different patches have different phases. The energy of such metastable configurations is 
only sligthly larger than that of the ground state. This scenario is clearly reminiscent of 
glassy systems. 



II. THE SYSTEM 



The BH model is currently realized in laboratories in terms of ultracold bosonic gases 
trapped op.ca. lattices Qy y Hy Q. fl- It. Ha.nton,a„ .eads 
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where M is the lattice size, a^, a| and rii destroys, creates and counts bosons at lattice 
site i, respectively. The local boson-boson interaction Ui, on-site potential Vi and hopping 
amplitude across neighbouring sites Ji^h can be tuned by varying controllable experimental 
parameters such as the atomic scattering length (via Feshbach resonance) and the strength 
and setup of the electromagnetic fields giving rise to the optical lattice. Disorder can be 
in principle introduced in all these parameters. While the case with random Vi has been 
mostly studied, several works have addressed models with random Ui or Ji^h and the phase 



diagrams thereof 
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2l|. For the sake of simplicity, in the following we will assume 



the local (repulsive) interactions and the hopping amplitudes to be constant throughout the 
lattice, Ui = U = 1, Ji^h = J > (that is, the on-site repulsive strength is our energy scale). 
Also, being interested in the phases of the system, we will consider a situation as closer to 
the thermodynamic limit as possible. Hence we will assume periodic boundary conditions, 
and that the harmonic confining potential typical of experimental system is so weak to be 
safely ignored. In summary, the only site dependent quantity in our system will be the 



local potential vt. We will consider two realization of disorder: local potentials uniformly 
distributed in [—A, A] and binary disorder. In the latter case Vi assumes the values A and 
with probabilities p and 1 — p. 



A. Zero-temperature phase diagram 

Let us now briefly recall the main features of the phase diagram of Hamiltonian ([T]). The 
different phases are characterized by the properties of the ground-state |\E') or of the low-lying 
sector of the spectrum. Standard quantities used to characterize the phase diagram are the 
energy gap between the ground and the first-excited state, the compressibility k = M^^S^A/", 
where A/" is the total boson population controlled by chemical potential /i, the condensate 
fraction /c, i. e. the maximal eigenvalue [^^j of one-body density matrix pj = ($|a|a/i|$), 
the SF fraction /s measuring the response of the system under the action of an infinitesimal 
to a phase twist or boost (see, e. g., {2^). 

The phase diagram is usually drawn in the fi/U — J/U plane, and it is mostly taken 
by the SF phase, characterized by a gapless spectrum, finite compressibility k > 0, finite 
condensate and SF fractions /s, /c > 0. The small- J/f/ region of the fi/U — J/U plane is 
occupied by a series of MI lobes, where the system is characterized by a gapped spectrum 
and vanishing k, /s and fc- In the presence of disorder further phases may appear [llj], the 
best-known being the BG. This was originally characterized as an insulating (non-SF) yet 
compressible (gapless) phase l|, l2|. Later it was shown that, at least at the mean-field level, 
the BG is marked by finite k and fc and vanishing /s 23|. 



B. Gutzwiller Mean-Field Approximation 

The Hilbert space of Hamiltonian ([1]) is infinite even on a finite system. However, one can 
take advantage of the total number conservation arising from the commutation of H and the 
total number operator = rij. This allows to address the properties of the BH Hamil- 
tonian separately in each fixed-number Hilbert subspace, whose size is d{J\f, M) = ("^^jjf 
for A/" bosons on a size-M lattice. Such a size makes exact diagonalization prohibitive also for 
relatively small lattices at fillings of the order of unity. Larger systems can be addressed by 
resorting to Complex and computationally demanding numeric simulations, such as quantum 
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Monte-Carlo, Density-Matrix Renormalization Group and Time-Evolving Block Decimation 
algorithms. 

An alternative approach giving satisfactory qualitative results at much lower computa- 
tional cost is provided by the site- decoupling (Gutzwiller) mean-field approximation. Ba- 
sically, it relies on the posit aja^ = aja^ -|- a*ak — <y*j<yk, where aj = (aj) and (■) denotes 
expectation value on the ground-state |\E') of the system. As a result, the mean-field BH 
Hamiltonian becomes the sum of "decoupled" on-site terms, H = J^j'^jy 



U = ^^ji^j - 1) + (^'j - /^) - J {(^hi + <^ilf) (2) 



where 7^ = X],- ^jhO^h- This in turn results into the mean-field ground-state being a product 



■'3 

of on-site states 



M 



j=l n=0 

where \Q) is the vacuum state, aj\Q) = 0. Note that the constraint on the 7j's preserves 
some degree of coupling among neighbouring sites. Also, we recall that finding the ground- 
state of Hamiltonian Ti is equivalent to finding the ground-state (lowest-energy fixed-point) 
of a set of coupled nonlinear dynamic equations more general than the discrete Gross- 



Pitaevskii Equations 28|, |29|. Such ground-state is fully determined by the set of on-site 



order parameters {aj} such that {iljj\aj\ipj) = aj, where ipj is the ground-state of the on-site 
mean-field Hamiltonian ([2]). Quantities fc and /s, as well as the other quantities generally 
employed in the characterization of the different phases of the model, are straightforwardly 
determined from once set {aj} is known (see, e.g., 24|). 

We recall that, unlike Eq. ([1]), the mean- field Hamiltonian does not commute with the 
total number of bosons: [H, N] 7^ 0. Hence the chemical potential fi appearing in the local 
mean-field Hamiltonians ([2]) has to be properly set in order to attain the desired boson 
population Af = '^j{'ipj\nj\ipj). 



III. RESULTS 



In the following we address a BH model on a 50 x 50 lattice with periodic boundary 
conditions and Uj = U > 0, Jjh = J > 0. As to the on-site potentials vj, we consider 
two different realizations of disorder, namely random potentials uniformly distributed in the 
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interval [—A, A], with A = 0.2 and potentials assuming the value A = 0.5 with probability 
p ^ 5 • 10~^. 

For both the realizations of disorder we set the Hamiltonian parameters so that the (final) 
set {ctj} displays a patchy pattern, where blotches with aj ^ are separated by regions 



where a^-'s virtually vanish. This has the effect of quenching the SF fraction |24j. l30l| . Con- 
versely, the condensate fraction is small yet significant, being basically /c = 
The resulting system is not SF, yet it retains some degree of coherence. Since this phase 
arises in the presence of disorder and is neither SF nor MI, we identify it with the BG 24 1. 
This is illustrated in Figs. [1] and [3l referring to the uniform and binary distribution of the 
disordered local potentials, respectively. 

Technically, the system is initialized by choosing a set \oLj\ which clearly does not satisfy 
the above described constraint. Subsequently, the sites of the system are addressed singularly 
and the relevant local order parameter is changed in order to minimize the energy while 
meeting the constraint. At each iteration all of the sites undergo one of these local moves. 
During this process the chemical potential is adjusted in order to achieve the desired boson 
population. We assume that the system has reached convergence when maxj (1^^^^ ~ '-^j I) — 
2.5 ■ 10~^, where s and s + 1 denote subsequent steps in the iteration process. When this 
condition is achieved, the constraint on the total population is met with a much higher 
precision. 

In principle the local order parameters are complex numbers, but we clearly observe 
that the state minimizing the system energy is real. More precisely, the lowest energy is 
attained by many configurations characterized by the same local moduli and a different 
global phase, = \aj\e^'^. Much interestingly, we also find a large number of very low lying 
states characterized by basically the same patchy \aj\ distribution as the ground-state(s) 
and a non-trivial phase distribution. The phase of the local parameters is not constant 
throughout the lattice, but only inside each non- vanishing \aj\ patch. Different patches 
feature in general a different "global" phase. 

Let us now discuss our results in more detail. As we say. Fig. [1] addresses the case of 
uniformly distributed local random potentials. The potential landscape Vj G [— 0.2f/, 0.2U] 
is shown in the upper left panel of the figure. The corresponding colorbar is the topmost 
one in the lower left quadrant of the figure. We choose to load the 50 x 50 lattice with 
A/" = 2650 bosons, corresponding to a filling slightly larger than unity. We find that at 
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J/U = Q - 10~^ the ground state is such that most of the lattice sites have a local population 
of one single boson, with a very small deviation. This is clear from the topmost left panel, 
showing the local population, where the colormap ranges from to 2 and unitary population 
corresponds to a green hue. The local order parameter aj virtually vanishes at these "highly 
squeezed" sites, as it is illustrated in the lower right panel of Fig. [H The corresponding 
colormap is the same as above, but the color range is aj G [0, 0.527]. The sites with virtually 
vanishing local order parameters have been assigned a black hue in order to highlight the 
blotchy pattern of the nonzero a^-'s. The patches of nonzero a^'s add up to a small yet 
defininitely nonvanishing condensate fraction, whereas the "sea" of vanishing local order 
parameter results in a virtually zero /s. As a side effect, this "sea" allows for the existence 
of very low-energy complex-a^ configurations characterized by the same local population 
and \aj \ pattern as the ground-state, but different phase pattern. Figure [2] shows such phase 
patterns for four of these low-lying configurations. The color key for these figures is encoded 
in the lowermost colorbar in the lower left quadrant of Fig. [H which ranges from to 27r. 
The plot in the same quadrant shows the relative energy difference between the ground-state 
(0) and the four configurations in Fig. [2] (1-4). Note that the islands of nonzero local order 
parameters in Fig. [2] have the same boundaries, and that all the sites belonging to the 
same island have virtually the same phase —i\og{aj/\aj\). Different islands in the same 
configuration or even the same island in different configurations have in general different 
"island-phases", recognizable as different colors. 

These low-lying states appear to be quite robust to our algorithm searching for the 
ground-state of the system, and in this sense we refer to them as metastable states. Our 
algorithm gets stuck at these low-energy configurations virtually indefinitely and, in this 
sense, we refer to them as metastable states. A similar behaviour has been recently observed 
in related systems, even in the absence of disorder [3l|, l32|. 

Figure [3] shows a similar situation on a 50 x 50 lattice where a binary disorder is present. 
In 128 out of the 2500 total sites the local potential equals 0.5U, whereas it is zero in the 
remaining sites of the lattice. The local potential landscape is illustrated in the top left 
quadrant of Fig. [31 It could describe the effect of 128 atoms of a second atomic species, 

lopping amplitude, which interact 



whose position is "frozen" due to a greatly quenched 
repulsively with the atoms of the bosonic species [lol . 
choose a filling close to unity, yet slightly lower, Af = 2455, and J/U = 1.8 • 10~^. Once 



20l |. For the latter we once again 



7 



again the local population of the ground state is 1 with a very small variance at most of 
the lattice sites. The remaining sites, which are those where the impurities are localized, 
and possibly their nearest neighbours, have a smaller population and a larger population 
variance. This is clear from the upper right panel in Fig. [3] (the relevant colormap is the 
same as in Fig. [H with the only difference that it ranges in [0, 1] so that in this case unitary 
filling corresponds to a dark red hue). The lower left panel of Fig. [3] shows that nonzero local 
order parameters appear at the impurity sites and possibly at their nearest neighbours. The 
remaining sites of the lattice have a virtually zero local order parameters, which produces 
the same effect as discussed above. Once again we find low-lying states characterized by 
virtyally the same blotchy distribution of the (modulus of the) local order pameter, but 
different phases on different nonvanishing-aj islands. This is clearly illustrated in Fig. HI 
whose color key is the same as Fig. [2j 

As J /U is increased the nonzero- local-order-parameter islands swell, and eventually melt 
together. This allows a superfiuid current to flow through the system, and the ground-state 
is not any more in a BG phase, but rather in a SF phase. Under these circumstances we find 
no low-lying states characterized by a random phase pattern of the local order parameter. 
This either means that such states do not exist or that they are not robust to our algorithm 
as those observed in the BG phase. Even starting our algorithm with local order parameters 
exhibiting a random phase pattern ends up in a "real" ground-state (modulo a global U(l) 
symmetry). When convergence is achieved the phases of the local order parameters are the 
same to a good approximation. 

IV. DISCUSSION 

We consider a two-dimensionl BH model with random local potentials, focusing on two 
different disorder distributions, uniform and binary. We show that in the BG phase many 
low-energy metastable configurations exist with a slightly larger energy than the ground- 
state. These configurations are characterized by the same blotchy pattern of local order 
parameters as the ground state. However, while in the latter the local order parameters 
have the same global phase, in the latter each nonvanishing-Oj island is characterized by a 
different constant phase. This situation is strongly reminiscent of the energy landscape of 
glassy systems. 
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We observe that the random phase pattern characterizing the low-lying states in the BG 
phase could account for the experimental observation that the BG appears to be lacking 
coherence 18|]. Indeed this is apparently in contrast with our characterization of the BG 



as a non-superfluid yet coherent phase. Clearly, one should take into account the fact that 
the spatial arrangement of the coherent blotches characterizing the BG is very complex, 
following from the random nature of the local potential pattern. The interference pattern 
produced by a large number of coherent yet randomly localized sources is very likely not so 
different from that of a set of incoherent sources. One could in principle think of "factoring 
out" the geometric contribution to the interference pattern. For instance, in the binary 
disorder case one could observe that the location of the coherent sources is basically the 
same as that of the random impurities. If the location of the latter is somehow known, it 
should be possible to recognize whether an interference pattern is actually incoherent or it 
looks so due to the random location of the coherent sources. 

However the incoherence of the interference pattern could arise from the "glassiness" of 
the BG. Due e. g. to imperfect cooling the system could be on a slightly excited state 
instead of the ground-state. If this is the case, correlating the interference pattern with the 
distribution of coherent source could prove useless, because anyway each of these sources 
could be characterized by a different unknown phase. 
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FIG. 1: Uniformly distributed random potential (A/[/ = 0.2) on a 50 x 50 lattice with periodic 
boundary conditions containing N = 2650 bosons with hopping amplitude J/U = 6 • 10"'^. Upper 
left panel: random local potentials; Upper right panel: (rij); Lower right panel: |(aj)p; The 
colormaps for these data are shown in the lower left quadrant of the figure. From top to bottom: 
local potentials, local density and order parameters (color ranges [—0.2U,0.2U] and [0,0.527], 
respectively), local phases, referring to Fig. [2] (color range [0, 27r]). The plot below the colorbars 
shows the relative energy differences between the real (0) and complex configurations (from 1 to 
4, displayed in fig. [2]). 
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FIG. 2: Local phase patterns for the same situation discussed in Fig. [H The four configurations 
have been obtained starting from different initial random phase patterns. The phases range from 
to 27r as described by the lowermost colorbar in Fig. [TJ Black areas correspond to region of 
vanishing |(aj)p. As discussed in the text, the {tij) and |(aj)p distributions of these complex 
configurations are virtually undistinguishable from those of the real configuration, displayed in the 
rightmost panels of Fig. [TJ 
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FIG. 3: Binary random potential (A = 0.5, 128 impurities) on a 50 x 50 lattice with periodic 
boundary conditions containing = 2455 bosons with hopping amplitude J = 1.8 • 10~^. The 
panels have the same meaning as in Fig. [TJ The color ranges for the local potentials (top colorbar), 
densities and order parameters (middle colorbar) are [0,0.5], [0, 1],[0, 0.263], respectively. 
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FIG. 4: Local phase patterns for the same situation discussed in Fig. [3l The four configurations 
have been obtained starting from different initial random phase patterns. The phases range from 
to 27r as described by the lowermost colorbar in Fig. [3j Black areas correspond to region of 
vanishing |(aj)p. The situation is analogous to that of the uniformly distributed random local 
potentials discussed in figEJ 



14 



